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1.  INTRODUCTION  AND  OVERVIEW 


N 

DELFIC  (DEfense  Land  £a!1out  £nterpretative  £ode)  is  intended  for  re¬ 
search  in  local  nuclear  fallout  prediction  and  to  serve  as  a  standard 
against  which  predictions  by  less  capable,  production-oriented  codes  can 
be  judged.  By  local  fallout  .We  means  the  intensely  radioactive  material 
which  falls  to  the  ground  within  several  to  several  hundred  miles  of  ground 
zero,  depending  on  the  size  of  the  explosion.  The  code  is  essentially  open- 
ended  with  regard  to  input  data,  it  is  highly  flexible  in  that  it  offers 
many  options  that  would  not  be  available  in  a  production-oriented  code,  and 
it  strives  to  include  as  much  of  the  physics  of  fallout  transport  and  ac¬ 
tivity  calculation,  without  resorting  to  short  cuts,  as  is  practicable. 

Calculation  begins  at  about  the  time  the  fireball  reaches  pressure 
equilibrium  with  the  atmosphere.  Rise,  growth  and  stabilization  of  the 
nuclear  cloud  is  computed  by  a  dynamic  model  that  treats  the  cloud  as  an 
entraining  bubble  of  hot  air  loaded  with  water  and  contaminated  ground 
material.  The  fallout  particle  cloud,  including  the  stem,  is  formed  during 
the  cloud  rise.  This  calculation  requires  specification  of  a  vertical  pro¬ 
file  of  atmospheric  data:  pressure,  temperature,  humidity  and  wind;  thus, 
the  height,  dimensions  and  vertically  sheared  horizontal  displacement  of  a 
particular  cloud  are  determined  by  the  atmospheric  stability  and  winds. 

After  cloud  stabilization,  representative  parcels  of  fallout  are  trans¬ 
ported  through  an  atmosphere  that  is  defined  by  input  data.  The  user  may 
specify  a  single  vertical  wind  profile  and  assume  a  steady  state.  He  may 
specify  any  number  of  wind  profile  updates.  He  may  resolve  the  transport 
space  in  the  horizontal  and  specify  multiple  wind  profiles  defined  at  dif¬ 
ferent  geographical  locations,  in  which  case  winds  in  the  cells  of  a  three- 
dimensional  space  grid  are  determined  by  an  interpolation  procedure. 
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During  transport,  fallout  parcels  are  expanded  in  the  horizontal  by 
ambient  turbulence.  Turbulence  data  may  be  input  along  with  the  winds, 
but  since  these  data  are  rarely  available,  they  can  be  calculated  by  the 
code. 

To  account  for  effects  of  vertical  wind  shear  on  the  dispersion  of 
individual  parcels,  tops  and  bases  of  the  parcels  are  transported  to  ground 
impaction  separately,  arid  then  recombined.  The  impacted  parcels  are  dis¬ 
tributed  over  the  ground  via  a  bivariate  Gaussian  function. 

Any  or  all  of  sixteen  unique  quantities  computed  from  the  deposited 
fallout  may  be  mapped.  Radioactivity  is  calculated  rigorously  for  any  time 
by  summing  exposure  or  exposure  rate  contributions  from  all  nuclides  in  the 
mixture  of  fission  products.  Any  of  twelve  different  types  of  fission  may 
be  specified.  Induced  activity  in  soil  material  in  the  fallout  and  in  238U 
may  be  accounted  for. 

The  remainder  of  this  presentation  is  partitioned  into  four  parts,  in 
accord  with  the  code  structure.  We  begin  with  Initialization  and  Cloud  Rise, 
proceed  to  Atmospheric  Transport,  then  to  Activity  Calculation,  and  finally 
we  discuss  Map  Preparation.  The  last  chapter  presents  some  validation 
resul ts. 

Volume  II  of  this  set  is  a  user's  manual. 
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2.  INITIALIZATION  AND  CLOUD  RISE 


2.1  INITIALIZATION 

2.1.1  Time,  Temperature  and  Mass  of  Fallout 

The  cloud  rise  calculation  is  initialized  at  roughly  the  time  when 
the  fireball  reaches  pressure  equilibrium  with  the  atmosphere.  Time  and 
temperature  as  functions  of  yield  and  height  of  burst  were  determined 
from  measurements  made  from  cinefilms  of  nuclear  cloud  rise  and  from  radia¬ 
tion  measurements  of  fireballs.  Details  are  given  in  ref.  2.  Initial  time, 
t^  (seconds)  is 


‘i  -  56  ,;2mr°'30 


(2.1.1) 


where  the  time  of  the  second  temperature  maximum,  t2m  (seconds),  is 

t2m  =  0.037  (l.216X/18°)w^0,49  "  °-07*/l80)j  180,  (2.1.2) 


W  is  yield  (KT)  and  X  scaled  height  of  burst  (ft  Kl"1/3)  (Subroutine  TIMEE). 
Initial  gas  temperature,  T..  (°K),  is 


T.  =  Kl 


fc) 


+  1500 


(2.1.3) 


in  which 


K  .  5980  (1.145)i/,8°  039 5+0. 0264X/ leo) 


and 


n  -■=  -0.4473W 


0.0436 
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Temperature  cf  condensed  phase  matter,  T.  _•  (°K),  is  conjectural  hut  important 
since  the  energy  required  to  heat  it  must  he  taken  from  that  available.  The 
following  specification  has  been  found  to  work  satisfactorily 

T  .  =  200  log10  (W)  +  1000  (2.1.4) 

where  W  is  in  kilotons.  (Subroutine  TEMP) 

Mass  of  fallout,  ms  (kg),  (i.e.,  soil  plus  weapon  debris)  in  the  cloud 
at  ti  is  given  for  a  subsurface  burst  with  scaled  depth  of  burst 
d(ft  KT"1/3-4)  by 


ms  =  kdW3/ 3- 4  R(d)2  D(d),  0  <  d  <  20 


(2.1.5) 


where 


R(d)  =  112.5  +  0. 755d  -  9.6  x  10"6d3  -  9.11  x  10_12d5 
D(d)  =  32.7  +  0. 85 Id  -  2.52  x  I0“sd3  +  1.78  x  10"10d5, 


and  for  surface  and  above  surface  bursts  with  scaled  height  of  burst  A(ft  KT"1/3*4) 
by 


ms  =  kAW3/3-4(l80  -a)2  (360  +  A),  18J  i  A  s  0.  (2.1.6) 

The  constants  k^  and  k>A  were  determined  from  an  analysis  of  Teapot  Ess  fallout 
to  be  2.182  and  0 . 07741 ( kg  ft3)  respectively.  The  scaling  equation  for  sub¬ 
surface  bursts,  eq.  (2.1.5),  is  based  on  Nordyke's  scaling  function  for  high 
explosive  cratering  results,14  and  the  above-surface  scaling  function,  eq. 
(2.1.6),  is  based  on  the  volume  of  intersection  of  a  fireball  with  the  ground.2 
(Subroutine  MASS) 
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When  A  >  180  (ft  KT"1/3*4),  we  assume  a  "pure  airburst"  (i.,e.,  the  fire¬ 
ball  does  not  touch  the  ground),  and  usually  no  local  fallout  is  produced. 

In  this  case  we  take 

t2  =  0.045W0-4*  (2. 1.2') 

and  t.j  is  then  calculated  from  eq.  (2.1.1).  Temperature  is 

Ti  -  K' (t./t2m)n  +  1500  (2.1.3* ) 

where 

K*  =  6847W-0*0131 

and  n  is  as  for  eq.  (2.1.3).  Mass  of  weapon  debris  is  taken  to  be  90.7  kg.* 
(Subroutine  AIRBRS) 


2.1.2  Soil  Sol idification  Temperature 

Of  great  importance  to  the  activity  calculations  (sec.  4.2.2)  is  the  time 
after  detonation  when  soil  particles,  which  will  constitute  the  bulk  of  the 
local  fallout,  cease  to  absorb  radioactive  fission  products  on  their  surfaces. 
In  DELFIC  this  is  taken  to  be  the  time  the  average  cloud  temperature  reaches 
the  so-called  soil  solidification  temperature.  The  user  may  specify  this 


Though  during  the  cloud  rise  calculation  all  soil  is  taken  to  be  in  a 
condensed  phase,  an  estimate  of  the  partition  of  the  soil  burden  between  vapor 
and  condensed  phases  is  made  and  printed.  Mass  of  vaporized  soil  at  t.,  m  . 
(kg),  is  computed  via  the  conjectural  equation  ’ 

Vi  ■  0.0001Sms(T1  -  Tm). 

where  T  is  3000  °K  for  siliceous  and  3100  °K  for  calcareous  soils  as  these 
are  defined  in  sec.  2.1.2.  (Subroutine  VAPOR) 


temperature,  or  on  default  of  input,  the  code  will  specify  it.  For  this 
purpose  we  distinguish  two  types  of  soil:  siliceous  soil  of  continental 
(e.g.,  Nevada  Test  Site)  origin  (2200  °K)S  and  calcarious  soil  of  coral 
(e.g.,  Pacific  Test  Site)  origin  (2800  °K).  On  default  of  soil  type  speci¬ 
fication,  the  code  selects  the  continental  type. 


2.1.3  Mass  and  Geometry  of  the  Cloud 

On  the  basis  of  considerable  experience  with  the  cloud  rise  model  we 
take  the  fraction  of  explosion  energy  used  to  heat  air,  soil  and  water  to 
their  initial  temperatures  to  be  45%  of  the  joule  equivalent  of  the  total 
yield,  W.  That  is,  the  energy  available  to  heat  the  cloud  contents,  H 
(joules),  is 

H  =  0.45  (4.18  x  1012W) 


where  W  is  in  kilotons.  To  allow  for  situations  where  substantial  amounts  of 
water  are  taken  into  the  fireball,  the  user  may  specify  a  number  cf>  which  is 
the  fraction  of  H  used  to  heat  air  and  soil;  the  fraction  (1  -  <f>)  then  is 
used  to  heat  condensed  water. 

If  the  amHent  temperature  at  the  initial  height  of  the  cloud  (eq. 

(2. 1.13))  is  T  .,  then  the  masses  of  air  and  water  in  the  cloud  at  t_. ,  m,  . 

G  j  I  i  a  j  I 

and  mw>-j(kg)  >  are 


♦ 


cs(T)dT 


cpa(T)dT+xe 


(T)dT 


(2.1.7) 
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(1  -  ♦)  H  -  m 


cs(T)dT 


i  * 

J  <yT)dT  +  L 


+  xA,i 


(2.1.8) 


The  total  initial  mass,  m. ,  <  id  volume,  V.,  of  the  cloud  are 


m.  =  m  .  +  m  .  +  m 

i  a,i  w,i  s 


(2.1.9) 


vi  ’  <ma.i  +n,w.i>  RJf/P  ' 


(2.1.10) 


Here  cs(T),  cpa(T)  and  cpw(T)  are  specific  heats  of  soil,  air  and  water  given 
by  (joules  kg-1  °K-1) 


cpa(T) 


=  946.6  +  0.19710T,  T  <  2300  °K 


cpa^T^ 


=  -3587.5  +  2.125T,  T  >  2300  °K 


cpw(T) 


=  1697.66  +  1.144174T 


cs(T) 


781.6  +  0.5612T  -  1.881  x  107/T2,  T  <  848  °K 


cs(T) 


=  1003.8  +  0.1351T,  T  >  848  °K; 


L  is  the  latent  heat  of  vaporization  of  water  or  ice,  taken  to  be  2.5  x  106 
and  2.83  x  106  joules  kg"1  respectively;  ^s  the  ideal  gas  law  constant  for 
air  (287  joules  kg"1  V1); 


xe  =  Vws^P) 


(2.1.11) 


13 


where  HR  is  relative  humidity  (fractional),  %  =  29/18  is  the  ratio  of  molecular 
weights  of  air  and  water,  P  is  ambient  pressure,  Pws  is  saturation  water  vapor 
pressure  (Pa)  given  by 

Pw$  =  611(273/Te)5- 13  exp[25(Te  -  273)/tJ  (2.1.12) 

where  Te  is  ambient  temperature  (°K);  and  T..  is  virtual  temperature  (see  the 
glossary  in  Appendix  D)  in  the  cloud  at  t.. 

The  initial  cloud  is  assigned  the  shape  of  an  oblate  spheroid  with 
eccentricity  0.75.  The  eccentricity  value  0.75  is  an  average,  with  standard 
deviation  0.08,  found  by  Norment  and  Woolf  for  ten  test  shots  that  cover  a 
yield  range  from  3.6  l(T  to  15  MT.15  Norment  and  Woolf  found  little  variation 
of  eccentricity  with  height,  particularly  for  small  and  medium  yield  shots,  up 
to  the  altitude  at  which  stabilization  and  final  horizontal  expansion  begins. 
For  a  cloud  with  eccentricity  0.75,  the  ratio  of  vertical  to  horizontal  cloud 
radii,  He/Re,  is  0.66144.  (Subroutine  CRMINT) 


2.1.4  Altitude  and  Speed 

Good  results  are  obtained  with  the  initial  cloud  center  altitude,  i 
(m  relative  to  sea  level),  given  by 


zi  =  zGZ  +  zHoB  +  90Wl/3  * 


(2.1.13) 


where  zQ^  and  zHqB  are  altitudes  of  ground  zero  and  height  of  burst  above 
ground  zero  and  W  is  yield  (KT). 

Initial  rise  speed,  u.(m  s"1),  is 

u.  =  1.2  SgR~  (2.1.14) 

I  t  )  I 

where  g  =  9.8  is  acceleration  of  gravity  and  R  .  is  initial  cloud  radius.  The 

C  j  1 
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form  of  this  equation  results  from  a  simple*  approximate  analysis  of  initial 
fireball  rise,  and  the  constant  1.2  is  chosen  to  fit  observed  data.  (Sub¬ 
routine  CRMINT) 


2.1.5  Atmospheric  Conditions 

Height,  radius  and  horizontal  location  relative  to  ground  zero  of 
the  stabilized  cloud*  are  sensitive  to  atmospheric  stability  and  the  ambient 
winds.  The  user  provides  single  vertical  profiles  of  atmospheric  temperature, 
pressure,  relative  humidity  and  wind  for  use  by  the  dynamic  cloud  rise.  They 
also  are  used  to  compute  particle  settling  and  advection  during  the  cloud 
rise  sucl,  that  by  stabilization  time,  the  particle  cloud  geometry  is  defined 
as  functions  of  particle  size  and  space  above  and  downwind  of  ground  zero. 
(Data  input  via  subroutines  ATMR  and  SHWIND ) 

Multiple  wind  profiles  which  may  be  used  later  for  atmospheric  transport 
(sec.  3.5)  are  not  used  here;  a  single  wind  profile  is  input  especially  for 
the  cloud  rise  calculation. 

In  addition  to  pressure,  temperature  and  humidity,  the  user  also  may 
input  air  density  and  viscosity,  but  in  lieu  of  their  input,  the  code  computes 
them  from  the  other  quantities.  Some  flexibility  of  input  is  allowed:  alti¬ 
tude,  temperature  and  relative  humidity  must  be  specified,  but  either  or  both 
of  pressure  or  density  may  be  specified. 


2.1.6  Particle  Size  Distribution 

The  user  may  specify  one  of  three  types  of  particle  size  distribution: 

1.  a  lognormal  distribution  of  number  of  particles  with  respect  to 
diameter,  specified  in  terms  of  median  diameter  and  geometric 
standard  deviation. 


Stabilization  occurs  .vhen  ambient  transport  and  dispersion  of  the  cloud  be¬ 
comes  dominant  over  internally  generated  rise  and  expansion. 
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2.  a  power  law  distribution  of  particle  mass  fraction  with  respect 
to  diameter,  specified  in  terms  of  the  power  law  exponent  and  a 
parameter  k/m  which  is  defined  below,  or 

3.  an  arbitrary  distribution  of  mass  fraction  with  respect  to  particle 
diameter,  specified  by  input  of  a  table  of  values. 

A  lognormal  distribution  is  selected  by  the  code  on  default  of  user  speci¬ 
fication. 


Lognormal  Distribution20 


The  lognormal  distribution  is  examined  in  some  detail  in  Appendix  A; 
here  we  present  a  brief  summary.  A  particle  distribution  is  said  to  be  log 
normal  if  the  number  of  particles,  dn(s),  in  the  diameter  range  ds  centered 
on  6  is  given  by 


(2.1.15) 


where  650  and  s  are  the  median  diameter  and  geometric  standard  deviation  of 
the  distribution. 

* 

For  DELFIC  calculations  the  user  may  specify  650  and  s  ,  or  on  default 
of  specification,  thp  code  assigns  the  values 

650  ~  0.407  ym 

s  =  4.0  , 

which  are  representative  of  Nevada  Test  Site  fallout  and  are  used  for  surface 
or  near  surface  bursts  (A  <  180  ft  KT"1/3*4),  while  it  assigns 

* 

See  Appendix  A  for  a  discussion  of  how  to  specify  <S50  and  s. 


650  =  0. 15  pm 
s  =  2.0 

for  pure  air  bursts  (A  >  180  ft  KT"1/3*4). 21 

A  unique  characteristic  of  the  lognormal  distribution  is  that  the  fre¬ 
quency  functions  for  particle  surface  area  and  volume*  are  simply  related  to 
the  frequency  function  for  number.  The  surface  area  frequency  with  respect 
to  diameter  is  obtained  by  replacing  £.ns50  in  eq.  (2.1.15)  with  i.n650  +  2(ans)2, 
and  the  volume  frequency  by  replacing  £,n65Q  in  eq.  (2.1.15)  with  Jtn6 50  + 

3(s,ns)2.  Thus  the  median  diameters  for  the  surface  and  volume  distributions 
corresponding  to  650  =  0.407  pm,  s  =  4.0  are  19  and  130  pm  respectively, 
while  those  for  650  *  0.15  pm,  s  =  2.0  are  0.39  and  0.63  pm.  (Subroutines 
I CM  and  DSTBN) 

Power  Law  Distribution 

Power  law  distributions  are  mathematically  meaningless  since  distribution 
functions  cannot  be  defined  for  them.  This  is  because  the  power  law  function 
is  not  properly  bounded  for  zero  argument.  Freiling  has  shown  that  fallout 
particle  distributions  that  have  been  represented  by  power  law  functions  can 
equally  well  be  fitted  by  lognormal  distribution  functions.22  The  implica¬ 
tion  of  Freiling' s  work  is  that  power  law  distributions  would  be  more  ac¬ 
curately  described  as  truncated  lognormal  distributions.  Nevertheless,  power 
law  distributions  may  be  useful  in  fallout  work. 

Define  the  power  law  frequency  as 

dn(s | k , X)  =  k6~Xd6  ,  (2.1.16) 

where  dm (5 j k ,X)  is  the  number  of  particles  in  the  diameter  range  d<S  centered 
on  5.  If  we  assume  spherical  particles  with  constant  density,  p,  we  have 


Note  that  volume  and  mass  distributions  are  equivalent  for  spherical  particles 
whose  densities  are  constant  over  the  distribution. 
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dfi 


(2.1.17) 


where  dF 


.<)  ■  * 
te  •>) 


3-X 


6m, 


is  the  fraction  of  the  total  fallout  mass,  ms>in  the 


diameter  range  <S  to  <5  +  d<$. 

The  mass  fraction  of  particles  in  the  macro  range  from  5.  to  <5.  is  ob- 

'  J 

tained  by  integration  of  eq.  (2.1.17)  between  these  limits  to  give 


aF 


0  <  X  <  4. 


(2.1.18) 


In  DELFIC  the  power  law  distribution  is  specified  by  the  exponential 
quantity  X  and  the  ratio  k/m.  Here  we  have  dropped  the  subscript  s  from  the 
mass  since  the  ratio  k/m  would  have  been  determined  from  one  or  more  samples 
of  fallout  rather  than  from  the  entire  mass  of  fallout.  Details  on  data 
analysis  and  a  description  of  the  distribution  in  histogram  form  are  given  in 
Appendix  B.  (Subroutine  DISTBN) 

Tabular  Distribution 

Mass  fraction  and  boundary  particle  diameters  for  each  particle  size 
class  of  a  distribution  histogram  are  specified  by  the  user.  The  central 
particle  diameter  for  the  size  class  is  taken  to  be  the  geometric  mean  of  the 
boundary  diameters.  (Subroutine  DSTBN) 

Size-Activity  Distribution 

The  user  may  specify  any  of  the  above  to  be  a  particle  diameter-activity 
fraction  distribution.  For  such  a  case  the  code  se^cts  an  activity  K  factor 
(Roentgens  m2  hr_]  KT-1)  to  match  a  user  specified  fission  type,  and  a  con¬ 
ventional,  rather  than  rigorous,  activity  calculation  ensues. 
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In  effect  this  procedure  also  is  used  for  pure  air  bursts,  since  DELFIC 
assumes  that  activity  is  uniformly  distributed  through  the  fallout  particles, 
regardless  of  their  size,  for  pure  air  bursts 


2.2  CLOUD  RISE 

2.2.1  Background 

The  DELFIC  cloud  rise  model  is  based  on  the  work  of  Huebsch16*17  as 
modified  by  Norment. 8 »13.  It  provides  a  dynamic,  one-dimensional,  entrain¬ 
ing  bubble  model  of  nuclear  cloud  rise,  which  is  based  on  a  set  of  coupled 
ordinary  differential  equations  that  represent  conservation  of  momentum,  mass, 
heat  and  turbulent  kinetic  energy.  The  nuclear  cloud  is  defined  in  terms  of: 
vertical  coordinate  of  its  center  (the  cloud  is  in  some  respects  treated  as  a 
point),*  cloud  volume,  average  temperature,  average  turbulent  eneryy  density, 
and  the  masses  of  its  constituents:  air,  soil  and  weapon  debris,  water  vapor 
and  condensed  water.  Cloud  properties  and  contents  are  taken  to  be  uniform 
throughout  the  cloud  volume.  The  differential  equations  are  solved  by  a 
fourth-order  Runge-Kutta  algorithm.  (Subroutine  RKGILL)  Complete  simula¬ 
tions  use  very  little  computer  time. 

Norment  presents  results  of  a  validation  study  of  the  model.13 

2.2.2  Differential  Equations  (Symbols  are  defined  in  the  Glossary  found  in 
Appendix  D. ) 


Momentum 


e‘-  1  g  - 


I_  o*  +  I 

T*  6  +  m  dt  u 


(2.2.1) 


Effects  of  this  model  limitation  on  rise  simulation  of  very  large  clouds  are 
discussed  in  sec.  4  of  ref.  13. 


Terms  (a),  (b)  and  (c)  represent  forces  due  to  buoyancy,  eddy-viscous  drag 
and  entrainment  drag  respectively. 

k  is  a  dimensionless  parameter  that  has  been  calibrated  to  give  satis- 

2 

factory  simulation  results  as: 

k^  =  max  [a. 004,  min(0.1,  O.IW'1/3)]  .  (2.2.1a) 

In  long  form,  this  equation  means: 

k2  =  0.1  ,  W  <  1  KT 

k2  =  O.IW-1/3,  1  <  W  <  15,625  KT 

k2  »  0.004  ,  15,625  <  W  KT  . 

The  characteristic  velocity,  v,  is 

v  =  max(u,  \/2E)  .  (2.2.1b) 


»'.I»||  I  m0h 


where  c  (T)  ■  B'c  (T)  +  (1  -  B')cs(T)  and  where  cp(T)  =  (cpa(T)  + 
xcpW(T)]/(l  +  x^‘  Terms  (a)*  (b)  anc*  (°)  account  for  the  effects  of 
adaibatic  expansion,  entrainment  and  dissipation  of  turbulent  energy  to 
heat  respectively. 

Wet  (water  is  condensed) 


,  L2xe 

T  - 'V»» 


(T  -  Te)  + 


L (x  -  xeM  1 


mg'  dt 


Te  P  V  a  /  P 


(2.2.3W) 


Turbulent  Kinetic  Energy  Density 


2k2  I*  0* 

Te 


(b) 

u2  1  dm 
2  m  dt 


(2.2.4) 


Terms  (a)  and  (b)  account  for  turbulent  energy  generated  by  the  mean  flow  via 
eddy-viscous  drag  and  momentum-conserving  inelastic-collision  entrainment 
respectively.  Term  (c)  represents  entrainment  dilution,  and  (d)  represents 
dissipation  to  heat,  where  the  turbulence  dissipation  rate,  e,  is 


k3(2E)3/2 


(2.2.4a) 


k2  and  v  are  given  by  eqs.  (2.2.1a,  b)  and  k3  is  a  dimensionless  constant 
of  value  0. 175. 
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Mass 


Dry  (water  is  not  condensed) 


(2.2. 5D) 


This  equation  is  based  on  observed  volumetric  growth  of  nuclear  clouds,  and 
term  (a)  contains  this  information.  The  other  terms  are  required  to  convert 
volumetric  to  mass  growth  rate.  Terms  (b)  and  (c)  account  for  temperature 
effects  caused  by  entrainment  and  conversion  of  turbulence  energy  to  heat; 
term  (d)  accounts  for  hydrostatic  expansion  of  the  rising  cloud. 

Wet  (water  is  condensed) 
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Equations  (2.2.5D  and  W)  account  for  mass  change  by  entrainment  of 
ambient  air  only.  The  net  change  must  also  account  for  fallout  at  rate  F 
to  give 


dm  dm |  F 

3‘  '  dtLt  "  ' 

(2.2.5) 

Good  results  are  obtained  when  the  dimensionless 

entrainment  parameter 

y  is  given  by 

* 

y  =  max  £max(0. 12, 

O.IW0-1)  ,  o.oiw1/3] 

(2.2.5a) 

In  long  form  this  equation  means: 

y  =  0.12 

W  <  6.192  KT 

y  =  O.IW0-1 

6.192  <  W  <  19,307  KT 

y  =  0.01W1/3  , 

19,307  <  W  KT 

• 

Soil  Mass  Mixing  Ratio 

(a) 

ds  1  1  +  x  1^  dm  I 

dt  g1  1  +  x  s  m  d t I 


ent 


(b) 

1  +  X  +  S  +  W  /  S  \  r 

m  s  + 


(2.2.6) 


Term  (a)  accounts  for  entrainment  dilusion,  and  term  (b)  accounts  for  fallout 
from  the  cloud. 

Water  Vapor  Mixing  Ratio 

Dry  (water  is  not  condensed) 


dx 

dt 


1  +  x  +  s  /  \  1  dm  I 

1  +  x„  'x  e  m  dt 


ent 


(2.2.7D) 
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Wet  (water  is  condensed,  x  is  the  saturation  mixing  ratio) 


-  =  (1  +  x/5)  Q  +  (1  +  x/O  -3-y  u  . 


at 


* 

V, 


(2.2.7W) 


Condensed  Water  Mixing  Ratio 


(a) 


dm 

dt 


ent 


(b) 

dx 

dt 


(2.2.8) 


Term  (a)  accounts  for  entrainment  dilusion  and  term  (b)  accounts  for  condensa¬ 
tion  of  water  vapor  in  the  cloud.  Fallout  of  condensed  water  is  neglected. 
(Subroutine  DERIV) 


2.2.3  Particle  Settling 

Particle  settling  speeds  are  needed  to  compute  fallout  from  the  cloud 
during  rise  and  later  to  define  in  detail  the  stabilized  particle  cloud  and 
to  compute  trajectories  of  individual  particles  from  the  cloud  to  ground  im¬ 
paction.  For  these  purposes  DELFIC  computes  still-air,  terminal,  gravity 
settling  speeds  of  spherical  particles.  Atmospheric  properties  at  the  particl 
altitude  are  obtained  by  interpolating  in  the  tables  of  data  supplied  by  the 
user  (sec.  2.1.5).  Particle  density  is  either  specified  by  the  user  or  set  by 
the  code  to  2600  kg  m-3. 

The  calculations  are  generalized  in  terms  of  the  dimensionless  Davies 
number,  NQ, 

N0  =  4p  (Pp  -  P)g  63/(3n2)  (2.2.9) 


where  p  and  pp  are  air  and  particle  densities,  6  is  particle  diameter,  n  is 
air  viscosity  and  g  is  acceleration  of  gravity.  Still-air,  terminal  settling 
speed,  f,  is,  according  to  Stokes'  law  for  small  Nq, 

f  =  r)NDs/ (24p6 )  ,  Nd  s  0.3261  .  (2.2.10a) 

For  intermediate  ND  use  Beard's18  equation 

f  =  Blexp  3.  18657  +  0.992696Y  -  0.153193  x  10“2Y2  -  0.987059  x  10“3Y3 

-  0.578878  x  10'3Y4  +  0.855176  x  10'4Y5  -  0.327815  x  10‘5Y6 

0.3261  <  Nd  <  84.175  (2.2.10b) 


where  Y  =  fl,n(Ng).  For  large  NQ  use  Davies'19  equations: 

f  =  nNDs  (4.166667  x  10“2  -  2.3363  x  10“4ND 

+  2.0154  x  10'6N2  -  6.9105  x  10"9N3)/(p6) ,  84.175  <  NQ  <  140, 

(2.2.10c) 

1  ogio  (fp«5/n)  =  -  1.29536  +  0.986X 

-  0. 046677X2  +  1.1235  x  10“3X3,  140  s  Np  <  4.5  x  107  ,  (2. 2. lOd ) 

where  X  =  logio(ND),  and  the  slip  factor  in  eqs.  (2.2.10a),  (2.2.10b) 
and  (2.2.10c)  is 

s  =  1.0  +  54.088nT1/2/(6P)  ,  (2.2.11) 

where  P  and  T  are  air  pressure  and  temperature.  (Subroutine  SETTLE) 
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2.2.4  Cloud  Shape  and  Volume 

Initially  the  cloud  is  given  the  shape  of  an  oblate  spheroid  of  ec¬ 
centricity  0.75  (sec.  2.1.3).  This  shape  is  maintained  until  the  cloud 
stops  rising,  after  which  the  vertical  thickness  of  the  cloud  remains  con¬ 
stant,  and  further  volume  growth,  if  any,  is  accommodated  by  lateral  expan¬ 
sion.  The  cloud  volume  at  any  time  is  computed  via 

V  =  RjVm/P.  (2.2.12) 

a 

(Subroutine  RKGILL) 


2.2.5  Effect  of  Wind  Shear 

To  account  for  effects  of  wind  shear,  simple  adjustments  are  made  to  the 
volume  terms  in  eqs.  (2.2.5D)  and  (2.2.5W).  Namely, 


S 

V 


/  S  ,  i  1.5 

v\7  ,H«r 


(2.2.13) 


where  vg  is  the  magnitude  of  the  wind  velocity  difference  between  the  top  and 
base  of  the  cloud  cap,  and  k6  is  a  dimensionless  constant  taken  to  be  unity. 
(Subroutine  DERIV) 


2.2.6  Termination  of  Cloud  Rise  and  Expansion 

There  are  two  normal  termination  switches,  and  stabilization  is  taken 
to  occur  as  soon  as  one  of  these  is  tripped: 


a.  The  R-RATE  switch  is  tripped  when 


R  y0. 014778 


and  u  =  0. 


(2.2.14a) 
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b.  The  U,  EK  switch  is  tripped  when 


E  <  max  [lO,  min  (23  +  9  1  og i qW ,  60)J  and  u  =  0.  (2.2.14b) 

In  long  form,  relation  (2.2.14b)  means  that  termination  occurs  when 


E  <  10 


W  <  0.0359  KT 


E  <  23  +  91og10W  ,  0.0359  <  W  <  12,915  KT 


E  <  60 


,  12,915  <  W  KT 


and 


u  =  0. 

Here  E  is  in  units  of  joules  kg"1  (i.e.  ,  m2s"2).  It  turns  out  that  most  low 
and  medium  yield  cases  are  terminated  by  this  switch,  but  megaton  yield  cases 
are  terminated  by  the  R-RATE  switch.  (Subroutine  CXPN) 


2.3  DEFINITION  OF  THE  PARTICLE  CLOUD 

The  set  of  coupled  differential  equations  described  in  sec.  2.2  are  inte¬ 
grated  from  time  t..  to  cloud  stabilization  such  as  to  provide  a  complete 
history  of  the  cloud  rise  and  growth.  During  this  calculation,  cloud  proper¬ 
ties  are  tabulated  at  intervals  in  an  array  CX,  the  cloud  history  table,  which 
appears  in  the  printed  output  (II,  p.  42).  (Subroutine  CXPN) 

To  define  the  particle  cloud,  the  cloud  rise  and  growth  is  resimulated 
for  each  particle  size  by  passing  through  the  CX  array.  Thus,  a  separate 
particle  cloud  is  defined  for  each  size.  At  time  t. ,  each  cloud  is  taken  to 
mo  /e  uniform  composition  and  a  cylindrical  shape  which  is  subdivided  in  the 
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vertical  into  KDI  cylindrical  parcels 
as  illustrated  in  Fig.  1.  The 
diameter  and  height  of  the  cylindri¬ 
cal  cloud  are  taken  equal  to  the  major 
and  minor  axes  of  the  spheroidal  cloud. 
KDI  is  either  specified  by  the  user  or 
assigned  the  value 


Figure  1.  Subdivision  of  an 
initial  cloud  cylinder  into 
four  parcels,  KDI  =  4. 


KDI  =  15  +  £n(W) 


(2.3.1) 


where  W  is  yield  in  kilotons. 

During  the  resimulated  cloud  rise  it  is  necessary  to  distinguish  between 
the  cloud  cap,  which  has  properties  as  recorded  in  the  C*  table,  and  the 
particle  cloud  which  is  being  calculated.  Gravity  settling  of  the  particles 
separates  them  from  the  cap  such  that  the  cap  and  particle  clouds  are  distinct 
at  all  times  subsequent  to  t.,  and  the  geometry  of  the  particle  cloud  as  dis¬ 
played  in  Fig.  1  applies  only  at  t.. 

As  the  CX  table  is  passed,  the  code  computes  the  height  and  diameter  of 
the  base  and  top  of  each  parcel  as  follows.  Particle  settling  speed,  f,  for 
actual  conditions  (in  or  below  the  cloud  cap)  is  calculated  (sec.  2.2.3),  and 
to  this  is  subtracted  an  upward  component,  u  ,  which  represents  the  influence 
of  the  cloud  cap  rise  on  the  particle. 


a.  In-cloud 


■  UB  +  (z  -  ZB> 


Gf^> 


(2.3.2a) 


b.  Below-cloud 


%  =  UB 


ZB  ~  z 
ZB  '  zGZ 


(2.3.2b) 


it 
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where  uT  and  ug  are  the  rates  of  rise  of  the  cloud  cap  top  and  base,  z  is 
the  parcel  base  or  top  altitude,  Zj  and  Zg  are  the  altitudes  of  cloud  cap 
top  and  base  and  zQZ  is  the  altitude  of  the  ground.  The  altitude  of  the 
parcel  base  or  top  is  computed  by  point-slope  integration  of  the  equation 


over  each  time  increment  of  the  CX  table. 

A  parcel  top  or  base  which  is  inside  the  cloud  cap  is  assigned  the  cap 
radius  as  recorded  in  the  CX  table.  When  the  parcel  base  or  top  falls 
through  the  base  of  the  cloud  cap,  it  is  assigned  the  cloud  cap  radius  at  the 
time  of  its  fallout,  and  it  keeps  this  radius  thereafter.  At  stabilization 
time,  the  base  and  top  of  each  parcel  is  examined  to  determine  if  either  or 
both  are  below  the  cloud  base.  If  not,  the  parcel  .is  assigned  the  cloud  cap 
stabilization  radius.  If  one  or  both  are  below  the  cloud  base,  the  parcel 
is  subdivided  into  n  parcels  where 


n  =  INT(Rt/Rb);  2  s  n  s  10  (2.3.3) 

and  Rt  and  Rg  are  the  radii  of  the  parcel  top  and  base.  Each  parcel  subdi¬ 
vision  is  taken  to  have  equal  vertical  thickness  and  equal  mass.  The  radius 
R  of  a  parcel  subdivision  at  altitude  z  is  computed  by  the  geometric  inter¬ 
polation  formula 

z  "  ZB  " 

ZT  "  ZB 

»  zb  -  z  -  ZT  (2.3.4) 

where  all  quantities  are  as  defined  above.  The  radius  of  the  ith  parcel  sub¬ 
division,  whose  top  and  base  are  at  altitudes  zi  and  z^_^ ,  is  computed  from 
eq.  (2.3.4)  for  the  altitude  of  its  center  of  mass,  zcm  ,  which  is 
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(Subroutine  RSXP) 

Each  parcel  of  fallout  also  may  be  subdivided  in  the  horizontal  as  con¬ 
trolled  by  input  parameter  I RAD.  While  the  code  retains  this  capability,  it 
has  not  proven  to  be  useful  and  will  not  be  discussed  further  here.  On  de¬ 
fault  of  input  of  IRAD,  no  horizontal  subdivision  occurs.  Details  are  given 
on  pp.  56-58  of  ref.  8. 

To  account  for  effects  of  advection  and  vertical  wind  shear  on  the  hori¬ 
zontal  disposition  of  the  particle  cloud,  the  CX  table  is  again  passed,  this 
time  backwards  for  each  fallout  parcel.  The  vertical  trajectory  of  the  parcel 
is  thus  recreated  and  its  horizontal  position  is  adjusted  according  to  the 
time  it  spends  in  each  wind  stratum  (sec.  2.1.5).  (Subroutine  WNDSET) 

Each  fallout  parcel  and  subparcel  is  defined  by  the  following  data 
which  are  passed  on  to  the  Diffusive  Transport  Module: 

1.  horizontal  space  coordinates  of  the  parcel  center 

2.  base  altitude 

3.  time  (stabilization  time) 

4.  vertical  thickness 

5.  radius 

6.  mass  (or  activity) 

7.  particle  diameter 

8.  volume 

The  aggregate  of  these  parcel  descriptions  represent  the  particle  cloud, 
which  includes  the  stem. 
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3.  ATMOSPHERIC  TRANSPORT 


3.1  GENERAL  DISCUSSION 

The  Initialization  and  Cloud  Rise  Module  passes  on  to  the  Diffusive 
Transport  Module  (DTM)  fallout  parcel  descriptions  as  discussed  at  the  end  of 
the  preceeding  section,  plus  vertical  profiles  of  atmospheric  pressure,  tem¬ 
perature,  humidity,  density  and  viscosity.  Wind  data  are  supplied  to  the  DTM 
as  either:  a  single  vertical  profile,  a  sequence  of  single  vertical  profiles 
used  to  update  the  windfield  and/or  multiple  profiles  to  account  for  varia¬ 
tion  of  wind  with  geographic  location.  If  the  latter  option  is  chosen,  the 
user  must  also  specify  a  three-dimensional  grid  which  spans  the  atmospheric 
transport  space  along  with  parameters  used  for  wind  field  interpolation,  and 
vertical  as  well  as  horizontal  wind  components  are  considered. 

During  transport  the  fallout  parcels  are  dispersed  by  ambient  atmos¬ 
pheric  turbulence.  A  scale-dependent  method,  corrected  for  effects  of 
gravity  settling,  is  used.  Turbulence  data  may  be  specified  by  the  user  in 
a  manner  similar  to  that  used  for  wind,  or  it  can  be  computed  by  the  code 
either  as  a  simple  function  of  altitude,  or  as  a  function  of  altitude  on  the 
basis  of  surface  data  specified  by  the  user. 

To  account  for  differences  in  wind,  turbulence  and  settling  speed  over 
the  vertical  span  of  individual  fallout  parcels,  the  tops  and  bases  of  the 
parcels  are  transported  separately  from  the  stabilized  cloud  to  the  ground. 
After  ground  impaction,  the  bases  and  tops  are  recombined  to  form  a  single 
deposit  increment  of  fallout  such  as  to  reflect  the  dispersion  between  them 
during  transport.  This  ensures  reasonably  accurate  fallout  prediction  even 
when  few  parcels  are  transported.  Furthermore,  n  +  1  transport  calculations 
are  required  for  n  parcels,  which  is  only  one  more  than  required  for  trans¬ 
port  of  parcel  center  points  alone. 


3.2  ADVECTION 

There  are  two  optional  modes  of  transport:  a  more  accurate  but  more 
time-consuming  "layer-by-layer"  method  which  computes  transport  through  each 
vertical  wind  field  stratum  in  a  stepwise  manner,  and  a  "quick"  method  which 
computes  transport  in  single  steps  from  initial  points  to  impact  (or  to 
time  or  space  grid  boundaries).  In  any  case,  particle  settling  speed  is 
computed  by  subroutine  SETTLE  for  each  vertical  stratum  using  atmospheric 
properties  for  that  stratum  for  each  particle  size.  As  noted  above,  fallout 
parcel  bases  and  tops  are  transported  independently. 


3.2.1  Layer-By-Layer  Transport 

For  a  general  case,  with  temporally  and  spatially  resolved  winds,  we  have 


J)LW.iki*Lto  -r 
IwUTKzTtT^TTzli 


(3.2.1) 


where  r(t  )  is  the  initial  horizontal  location  vector  of  the  parcel  top  or 
base  center,  z^  is  its  initial  altitude,  r(tj)  is  the  horizontal  location 
vector  of  the  parcel  base  or  top  after  impaction  on  the  ground  at  altitude  zg, 
fi(x,y,z,t)  is  the  horizontal  wind,  w(x,y,z,t)  is  the  vertical  wind  component, 
f(z)  is  particle  settling  speed,  and  the  summation  is  taken  over  vertical  wind 
strata  of  (variable)  depth  Az  through  which  the  particle  passes.  Horizontal 
spatial  resolution  of  the  wind  field  is  achieved  by  defining  wind  vectors  at 
the  centers  of  a  three-dimensional  space  grid.  The  code  is  provided  with  an 
interpolation  procedure  which  computes  the  grid  cell  wind  vectors  from  an 
arbitrary  set  of  spatially  resolved  input  data.  Temporal  resolution  is  pro¬ 
vided  by  updating  the  complete  wind  field  at  user  specified  intervals. 
(Subroutine  TRANP) 
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3.2.2  Quick  Transport 

Quick  transport  is  allowed  only  if  all  vertical  wind  components  are 
taken  to  be  zero.  Wind  data  are  presummed  over  the  vertical  wind  strata 
according  to 


j 

■  E  V2k  <3-2'2> 

k=l 


where  is  the  vector  sum  stored  for  the  jth  wind  stratum,  ^  is  the  wind 
for  the  kth  stratum  of  depth  Az^,  and  the  summation  is  taken  from  the  surface 
stratum,  k  =  1,  to  the  jth  stratum.  Similar  summations  are  stored  for  the 
turbulence  data  and  for  the  wind  vector  angle  relative  to  the  easterly 
direction.  (Subroutine  SUMDAT)* 

During  quick  transport,  we  require  the  average  wind  <  ft  >..  k  between  the 
bases  of  strata  j  and  k.  This  is  obtained  (Subroutine  GETDA)  via  equation 

<  U  >j ,k  =  ^k-1  ‘  ^j^~B,k  "  zB,j)  »  J  <  k’  (3.2.31 


where  zD  .  is  the  altitude  of  the  base  of  the  jth  stratum,  etc.  Similar 

D,J 

equations  apply  for  turbulence  and  angle.  By  making  use  of  a  precalculated 
table  of  deposition  times  from  each  stratum  base,  T  ,  we  have  for  a  simple 
case' that  does  not  involve  space  or  time  boundaries 


r(tj) 

U 


(t,  )  + 


Vn  • 


‘B.k,/f<k) 


X 


u  > 


1  ,k 


(3.2.4) 


where  f(k)  is  particle  settling  speed  and  is  wind  in  the  kth  stratum  and 
other  quantities  are  as  previously  defined.  For  cases  where  there  is  resol u 
tion  in  space  or  time,  partial  stratum  transport  terms,  such  as  the  middle 


ic 

The  wind,  turbulence  and  angle  arrays  are  replaced  by  these  sums.  When 
individual  quantities  are  desired,  for  example  for  layer-by-layer  transport, 
they  are  retrieved  by  calculation  from  the  summed  quantities.  (Subroutine 
GETDA) 
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term  on  the  right  side  of  eq.  (3.2.4),  are  computed  on  both  sides  of  a 
boundary  whenever  it  is  crossed.  Thus,  for  finely  resolved  cases,  this 
mode  may  lose  much  of  its  computational  advantage,  and  the  more  accurate 
layer-by-layer  method  should  be  used. 


3.3  TURBULENT  DISPERSION 

During  and  subsequent  to  atmospheric  transport,  fallout  parcels  are 
taken  to  have  Gaussian  distributions  in  the  horizontal.  As  the  parcel  is 
subjected  to  the  dispersive  action  of  the  turbulent  atmosphere,  the  vari¬ 
ance,  a2,  of  the  distribution  grows  accordingly.  The  pertinent  measure  of 
turbulence  is  turbulence  energy  density  dissipation  rate,  e(m2s"3),  which  is 
the  quantity  stored  in  the  turbulence  data  arrays.  In  accord  with  the 
Kolmogoroff- Batchelor  theory  of  scale  dependent  Lagrangian  diffusion  of  a 
puff  of  inert  material  by  homogenious  turbulence  as  formulated  by  Walton,23 
the  parcel  variance  at  time  t  is 

a(t)2  =  jo(ts)?-/3  +  ^ce1/3(t  -  ts)J  ;  a(t$)  <;  a(t)  <  (3.3.1a) 

o(t)2  =  o l  <2c(t  -  ts)(e/a|)l/3  +  2};  o(t)  >  a, 

(3.3.1b) 


where 


c 


[■  *(TTf 


(3.3.2) 


a2  is  the  parcel  variance  when  its  dispersion  rate  becomes  constant,  taken  to 
be  109m2  (this  value  is  quite  uncertain),  and  o(ts)  is  the  parcel  standard 
deviation  at  cloud  stabilization  time,  taken  to  be  one-half  the  parcel  radius 


★ 

Turbulence  data  arrays  are  congruent  with  the  wind  data  arrays. 
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supplied  by  the  Initialization  and  Cloud  Rise  Module,  c  is  a  correction  for 

particle  settling  derived  by  Csanady,24  where  <  f  >  is  altitude-averaged 

particle  settling  speed,  is  the  Lagrangian  time  scale  of  turbulence, 

is  Eulerian  length  scale  of  turbulence  and  n  has  values  1  and  2  for  downwind 

and  crosswind  turbulence  respectively.  Good  results  are  obtained  with 

t  =  n  * 

‘L  uE' 

Turbulence  data,  in  the  form  of  e  values,  may  be  specified  by  the  user 
in  the  same  manner  as  the  winds  are  specified.  However,  since  these  data 
rarely  will  be  available,  the  code  can  provide  representative  vertical  pro¬ 
files  of  horizontally  isotropic  e  via  the  method  of  Wilkins.26 

Wilkins  finds  that  e  can  be  approximated  by  a  simple  hyperbolic  function 
of  height  above  ground,  z,  by 


(ro/o)3'2 

‘  ■  ism 


(3.3.3) 


where  tq  is  atmosphe.  ic  surface  layer  shearing  stress,  p  is  air  density, 
k  =  0.35,  and  zq  is  surface  roughness  length. 


*Actually,  we  should  have  T^/D^  -  g/a  where  g  «  4  is  the  ratio  of  Lagrangian 
to  Eulerian  turbulence  time  scales  and  a,,  is  the  standard  deviation  of  verti- 
cal  turbulence.  According  to  Chen,25  for  example,  a •  =  5.39e1/3,  so  we  have 

c  =  [l  +  (0.74n  <  f  >  e-*/3)2]"^2  (3.3.2*) 


Since  e"1/3  is  typically  in  the  range  10  to  100,  this  correction  tends  to 
decrease  the  dispersion  rate  to  unreal istically  low  levels.  While  the  code 
provides  an  option  to  use  this  equation  (with  a  spatially  and  temporally 
averaged  e;  normally  we  take  =  1  as  noted  above. 
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Equation  (3.3.3)  also  can  be  written 


e 


k<2  +  v 


where  u*  is  surface  layer  friction  velocity  given  by 


(3.3.4) 


u*  =  k  us/  [  An(s/zo)  +  ipm]  . 


Here  u  is  surface  wind  speed  measured  at  height  s  (usually  s  =  10  meters), 
and  i|im  has  the  following  empirical  functional  forms27  of  height  and  Monin- 
Obukhov  Length,  L: 

a.  stable  conditions 

=  4.7  s/L  ;  L  >  0 

b.  neutral  conditions 

-  0  ;  L  =  ® 
rm 

c.  unstable  conditions 

^  "  *"[(?2  +  DU  +  l)2/8]  +  2  tan-1  (?)  -  f  ; 

5  =  (1  -  15  s/L)1/4;  L  <  0. 

The  user  can  input  to  DELFIC  values  for  us>  s,  zQ  and  1/L,  and  the  code  com¬ 
putes  e  via  eq.  (3.3.4)  at  the  altitude  of  each  vertical  wind  stratum.  Typical 
values  of  zq  are  given  by  Pasquill,28  and  1/L  values  art  given  by  Golder.29 
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For  many  predictions  it  is  sufficient  to  let  the  code  set  the  parameters 
in  eq.  (3.3.4).  In  this  case,  typical  (or  "average")  values  of  e  are  com¬ 
puted  via 

e  =  0.03/z  (m2s-3)  (3.3.5) 

3.4  DEPOSIT  INCREMENT  DESCRIPTION 

As  already  noted,  the  base  and  top  of  each  fallout  parcel  are  inde¬ 
pendently  transported  from  their  position  in  the  stabilized  cloud  to  ground 
impact.  Once  impacted,  the,.  -v  recombined  to  define  a  deposit  increment  of 
fallout  distributed  over  the  'impact  plane  by  a  bivariate  Gaussian  function. 

Figure  2  illustrates  how  the  distribution  parameters  of  the  base  (sub¬ 
script  b)  and  top  (subscript  t)  wafers  are  synthesized  to  yield  the  deposit 
increment  distribution  parameters.  The  subscript  quantities  „  and "7  aTc"* 
used  to  distinguish  alongwind  from  crosswind  standard  deviations.  The 
ellipse  orientation  angles,  and  o^,  are  the  wind  vector  angles  (relative 
to  east)  averaged  over  the  particle  trajectory  by  a  procedure  analogous  to 
that  used  for  the  wind  (eq.  (3.2.3)).  The  deposit  increment  ellipse  is  the 
large,  heavy  lined  ellipse.  The  deposit  increment  distribution  parameters 
are  computed  from 

a  =  a  r  c  t  a  n  (3.4.1)* 


* 

Since  the  FORTRAN  arctan  subroutine  yields  a  in  the  range  ±  ir/2,  it  is 
necessary  to  add  the  statement 

I F( XT  -  XB  •  LT  •  0.0)  ALPHA  =  ALPHA  -  S I GN ( 3  141592654 , ALPHA) 


to  extend  the  range  to  ±  ir. 


Figure  2.  Illustration  of  the  synthesis  of  a  deposit  increment  from 
the  standard  deviation  ellipses  of  the  top  and  base  of  a 
parcel  of  fallout. 
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A  a  +  and  A a  ,  are  computed  similarly,  and 

II  I  9  * 


XP  =  xb  +  (o„  ■  4o„,b>cos  “ 
yp  =  yb  +  <au  -  to  >b)s1n  c. 


(3.4.2) 


(3.4.3) 


The  time  recorded  with  the  deposit  increment  is  taken  to  be  the  arithmetic 
average  of  the  impact  times  of  the  base  and  top.  (Subroutine  ADVEC)  The 
center  coordinates,  xn,  y„,  standard  deviations,  a  ,  a  ,  orientation  angle, 

P  P  l  II 

a,  and  other  properties  received  from  the  Initialization  and  Cloud  Rise  Module, 
are  passed  along  to  the  Output  Processor  Module  where  contributions  in  the 
deposition  plane  are  computed  by  eq.  (5.2.1). 
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3.5  WINDFIELD  DESCRIPTION 


3.5.1  Single  Profile  Windfields 

Windfields  that  are  horizontally  homogenious  are  defined  by  single 
vertical  profiles  of  two-dimensional  wind  vectors.  (Vertical  components 
are  taken  to  be  zero.)  The  atmosphere  is  stratified  into  layers  based  on 
the  altitude  increments  at  which  the  data  are  specified  by  the  user. 
Limiting  horizontal  boundaries  of  the  transport  space  also  are  specified. 

Temporally  resolved  winds  may  be  of  this  type;  the  winds  are  replaced 
by  updates  at  specified  time  intervals.  The  only  restriction  is  that  the 
data  must  be  specified  at  the  same  altitudes  in  all  updates.  (Subroutines 
DATIN  and  ONEDIN) 


3.5.2  Multiple  Profile  Windfields 

Horizontally  nonhomogenious  windfields,  that  is  winds  that  vary  with 
geographical  location,  are  specified  by  multiple  vertical  profiles  of  winds 
which  may  be  three-dimensional.  (If  vertical  components  are  to  be  con¬ 
sidered  in  single  profile  windfields,  the  methods  of  this  section  must  be 
used.)  Multiple  profile  windfields  also  may  be  time  resolved. 

The  user  must  specify  a  three-dimensional  space  grid,  for  each  cell 
of  which  a  three-dimensional  wind  vector  is  to  be  defined  by  the  code  on 
the  basis  of  the  input  data.  A  table  of  vertical  strata  center  or  base  al¬ 
titudes  is  specified,  and  these  strata  levels  are  held  constant  throughout 
the  transport  space.  (Subroutine  LAVERS)  For  the  horizontal,  a  rectangular 
"control"  net  is  specified  in  terms  of:  the  coordinates  of  its  south-west 
corner,  a  (square)  mesh  spacing,  and  indices  which  represent  the  numbers  of 
mesh  units  to  the  east  and  north.  Each  control  net  mesh  may  be  quartei  ad, 
each  quarter  may  be  quartered,  and  so  on  (Vol.  II,  Fig.  A.l).  An  identical 
net  applies  at  each  vertical  stratum.  (Subroutine  GETUP) 
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Consider  a  general  wind  or  turbulence  component  whose  value,  D^, 
at  the  center  of  the  nth  atmosphere  space  cell  is  to  be  determined  from 
the  set  of  M  specified  components,  D^,  i  =  1,2,3,.,.,M.  Here  each  D.  is 
defined  at  an  arbitrary  point  in  the  atmosphere  space.  To  determine 
we  specify  an  integer  N  and  two  distances,  a  and  3,  and  then  the  code 
computes 


g.D. 
si  i 


(3.5.1) 


where  the  summation  is  over  the  N  data  with  the  largest  ,  where 

V  ♦  If/  y  ♦  X}  +  yy 
3i  '  H  71  Ia 

§(M) 


g2  -  xk  -  y£\ 

e +  +  Aj 


(3.5.2) 


Here  a  and  3  are  limiting  distances  in  the  vertical  and  horizontal  respec¬ 
tively.  Any  terms  in  the  g^  that  are  computed  to  be  negative  are  set  to 
zero.  This  weighting  method  is  similar  to  one  developed  for  similar  use  by 
Cressman. 30  (Subroutines  DATIN  and  TRIDIN) 
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4.  ACTIVITY  CALCULATION 


4.1  GENERAL  DISCUSSION 

Radioactivity  calculations  at  three  levels  of  sophistication  are 
available  in  DELFIC: 

1.  Rigorous  calculations  that  sum  contributions  from  all  nuclidas  in  a 
decay  chain.  Decay  from  detonation  to  the  requested  time  is  com¬ 
puted  for  each  nuclide  by  the  Bateman  equation  (eq.  (4.2.1)),  and 
each  activity  is  multiplied  by  a  unique  factor  to  convert  to  ex¬ 
posure  rate;  this  eliminates  the  need  to  use  a  power  law  decay 
f*:tor  which  is  applicable  to  the  mixed  fission  products.  Radio¬ 
chemical  fractionation  is  accounted  for  in  distributing  activity 
with  respect  to  particle  size,  and  activity  from  individual  mass 
chains  can  be  predicted  and  mapped.  A  conventional  K  factor  (i.e., 
exposure  rate  at  H  +  l  hour  over  fallout  from  one  kiloton  of  fission 
yield  spread  uniformly  over  unit  area)  is  not  used  by  this  method, 
but  one  can  be,  and  is,  calculated.  (Subroutines  PAM1,  PAM2, 

BATMAN,  GXPSR,  INDCD2 ,  MCHDEP  and  URAN) 

2.  The  rigorous  calculation  is  done  for  H  +  1  hour,  and  then  a  t"1*26 
decay  factor  is  used  for  subsequent  times.  This  option  is  offered 
to  save  computer  time  in  computing  integrated  exposure  (dose)  when 
time  of  arrival  of  individual  fallout  parcels  is  accounted  for.  A 
completely  rigorous  calculation  requires  repetition  of  the  entire 
procedure  for  each  fallout  parcel,  and  can  be  expected  to  use  a  lot 
of  computer  time, 

3.  For  pure  air  bursts,  and  particle  distributions  specified  in  terms 
of  size-activity  fraction,  K  factors  and  the  t"1*26  decay  factor 
are  used.  (Subroutines  PAM1A  and  PAM2A) 
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Twelve  types  of  fission  reactions  may  be  specified  for  rigorous  calcu¬ 
lation.  In  terms  of  the  input  parameter  FISSID  these  are: 


U233HE 

U235TH 

U235HE 

U233TH 

U235FI 

U233FI 

U238HE 

U238FI 

U238TN 

P239TH 

P239HE 

P239FI 

Here  U  and  P  represent  uranium  and  plutonium,  the  numerical  part  is  the  atomic 
mass  and  the  other  symbols  represent: 

HE  high  energy  neutron  fission 

FI  fission  spectrum  neutron  fission 

TN  thermonuclear  neutron  fission 

TH  thermal  neutron  fission 

For  the  non-rigorous,  level  3  calculation  K  factors  for  the  seven  most  common 
fission  types  in  the  left  column  are  available. 

Activity  calculations  are  controlled  by  the  Output  Processor  Module. 

4.2  RIGOROUS  ACTIVITY  CALCULATION 

The  following  presentation  is  substantially  the  same  as  section  2  of 
ref.  5;  it  outlines  the  operations  of  the  Particle  Activity  Submodule.  It  is 
given  in  terms  of  exposure  rate,  but  extension  to  exposure  should  be  self- 
evident.  Activity  from  individual  nuclides  also  can  be  computed  by  this 
method.  (Subroutine  MCHDEP) 
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4.2.1  Nuclide  Abundances 


Fission  Products 

The  fission  products  can  be  grouped  into  a  number  of  decay  chains, 
one  of  which  may  be  represented  by 


Ni  N2  +  N3  N 


m 


One  or  more  members  (including,  of  course,  Nj)  may  be  produced  directly  by 

-\jt 

the  fission  event,  and  each  member  except  Nm  decays  as  e  ,  while  each 

member  except  Nj  grows  by  e  .  The  complete  growth  and  decay  history 
for  each  member  of  a  chain  is  represented  by  the  Bateman  equation31 


n  n  -Xkt 

N„tt>  =LX>N°e 

i=l  k=l 


(4.2.1) 


where 


N  (t)  =  atoms  of  the  nth  member  of  the  chain  at  time  t, 

n 

N°  =  atoms  of  the  ith  number  of  the  chain  (i  <  n)  at  time  0, 


and  \ , 


=  disintegration  constant  of  the  kth  member  of  the  chain  (k  f  n), 


nCki 


k) 


(4.2.2) 


The  activity  in  disintegrations  per  second  is  given  by 


An(t)  =*nNn(t) 


(4.2.3) 
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The  number  of  disintegrations  that  occur  between  times  ti  and  t2  is  given  by 


ti t  t2 ) 


(4.2.4) 


which  is  obtained  by  integration  of  eqs.  (4.2.1)  and  (4.2.3).  The  Bateman 
equation  develops  a  singularity  whenever  This  difficulty  is  cir¬ 

cumvented  in  DELFIC  by  slightly  incrementing  one  member  of  such  a  pair. 
(Subroutine  BATMAN) 

Induced  Activity  in  Soil 

The  model  for  the  contribution  of  induced  activity  in  the  cloudborne 
soil  was  developed  by  Jones.32  Three  basic  assumptions  are  made: 

1.  All  neutrons  entering  the  soil  are  thermal ized  and  then  captured. 

2.  Only  those  neutrons  need  be  considered  that  are  seen  by  the  appar¬ 
ent  crater. 

3.  All  significant  soil  components  are  refractory  in  the  fractiona¬ 
tion  scheme. 

The  model  applies  only  to  Nevada  Test  Site  soil.  (Subroutine  INDCD2) 

Induced  Activity  in  Device  Materials 

The  model  currently  accounts  for  induced  activity  in  only  one  device 
component,  238U.  Neutron  capture  by  238U  produces  239U,  which  decays  to 
239Np  and  then  239Pu.  The  following  equations  are  used: 


1.  Activity  of  239U,  dis  s"1, 


Ai  =  A!Nie"Xlt 


(4.2.5) 
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number  of  neutron  captures 


N°  ■ 

Xi  =  disintegration  constant  of  239U,  s-1, 
and  t  =  time,  s. 

2.  Activity  of  239Np,  dis  s"1, 

XiX2N°  /  -X i t  -x2t\ 

A*  =  —xi  -  *  ;  • 

where 

X2  =  disintegration  constant  of  239Np,  s"1. 

3.  Disintegrations  of  239U  occurring  between  t*  and  t2, 


Disintegrations  of  239Np  occurring  between  ti  and  t2 


. . .  ir  " 


5.  Disintegrations  of  239U  from  ti  to  infinity. 


Mi 


=  Ni  e 


0  A i t i 


6.  Disintegrations  of  239Np  from  ti  to  infinity, 


M2 


Nl 


—  (\2e'Xltl  -  »!*"**)  ■ 


(4.2.9) 


(4.2.10) 


This  mass  chain  is  considered  completely  refractory  in  the  fractionation 
scheme.  (Subroutine  URAN) 


4.2.2  Activity  Versus  Particle  Size 

The  radial-distribution  model  of  Freiling33  categorizes  a  nuclide  as 
volatile  if  its  boiling  point  is  less  than  the  solidification  temperature  of 
the  soil  particles  (sec.  2.1.2).  The  nuclide  is  refractory  if  its  boiling 
point  is  greater  than  the  solidification  temperature.  Freiling  then  assumes 
that  a  volatile  nuclide  will  condense  on  the  particle  surfaces,  whereas  a 
refractory  nuclide  will  condense  uniformly  throughout  the  particle.  A  fission- 
product  mass  chain,  whose  composition  changes  with  time,  can  then  be  charac¬ 
terized  by  the  fraction  of  its  membership  that  represents  refractory  nuclides 
at  the  time  the  cloud  cools  to  the  soil-solidification  temperature.  This 
Freiling  rauo  is  referred  to  as  FR.  Freiling  has  postulated,  on  an  empiri¬ 
cal  basis,  that  the  specific  activity  of  fallout  particles  is  proportional  to 
the  (b.  -  1)  power  of  the  particle  size,  where  b.  is  /FD  for  mass  chain  i. 

1  1  K 


47 


Am. 


Frei ling34  has  treated  the  relationship  of  the  model  to  the  lop-normal 
distribution  of  particle  sizes.  In  DELFIC,  however,  the  model  is  general¬ 
ized  for  arbitrary  particle-size  distributions.  The  approach  can  be  formu¬ 
lated  as  follows: 

Let 


and 


6^  =  geometric  mean  diameter  of  kth  particle-size  class, 

=  number  of  particles  in  class  k, 

-Fj  =  total  number  of  equivalent  fissions  in  all  size  classes, 

Y.  =  fission  yield  of  ith  mass  chain, 

r  =  subscript  index  for  perfectly  refractory  mass  chain. 

The  equivalent  fissions  of  mass  chain  i  in  size  class  k  is 


W  ■  ViW 


(4.2.11) 


where 


W 


(4.2.12) 


For  the  perfectly  refractory  chain  (b..  =  1),  Eq.  (4.2.11)  takes  the  form 


MV  ■  ftVm<V  • 


where 


W  "  Nk5k  /  2  Nk6l  ’ 

i.e. ,  the  mass  fraction  in  the  kth  size  class. 
The  fractionation  ratio  is  given  by 


L  \  =  li  .  W 

\  1  -r/k  Yr  ¥¥ 


(4.2.13) 


(4.2.14) 


(4.2.15) 
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,  4-,  V-- 


(r^  A 


1 1 

Hvi 


(4.2.16) 


E  NkC 


To  obtain  an  algorithm  that  takes  advantage  of  the  fact  that  particle 
size  distribution  consists  of  a  table  of  fM(6^)  versus  6'k>  we  note  that 


n 


n 


Vk1  /SNkfik 

k=l  /  k=l 


i(wO")  • 

k=l  v  ' 


(4.2.17) 


Taking 


E.  *  1 


/i  (w  •  «bkrl)  • 


(4.2.18) 


we  can  rewrite  eq.  (4.2.16) 


(ri.r)k  =  %  4r'Ei  • 


(4.2.19) 


It  follows  from  eqs.  (4.2.11),  (4.2.15)  and  (4.2.19)  that 


W  ■  ft¥iV  f«<sk>  • 


(4.2.20) 


In  this  form  of  the  radial  distribution  model,  the  specific  activity  for  any 
mass  chain  that  is  not  perfectly  refractory  decreases  monotonically  with  in¬ 
creasing  particle  size.  Observations  of  fallout,  on  the  other  hand,  provide 
strong  evidence  that  the  specific  activity  tends  to  level  off  at  about  100  to 
200  uni.  To  account  for  this  effect,,  the  model  was  modified  in  the  following 
manner: 

The  approach  assumes  a  two-component  system.  One  component  obeys  the 

radial  distribution  model.  The  other  has  a  constant  specific  activity  over 

—  ( 

the  particle  size  distribution. 

Let 

R .  =  fraction  of  fissions  in  ith  mass  chain  that  obeys  radial 

distribution,  and 

Si  =  fraction  of  fissions  in  ith  mass  chain  that  appears  with 
constant  specific  activity. 

Then  Eq.  (4.2.20)  is  replaced  by 

W  ■  Vi  (RiEi4kr‘ +  si)  W  •  <4-2-21> 

Figure  3,  an  idealized  plot  from  fallout  observations,  shows  the  two  compo¬ 
nents.  (Subroutines  FRATIO  and  GXPSR) 


Figure  3.  Relationship  of  equivalent  fissions  to  particle  size 
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The  crossing  point  sets  the  criterion 


since 


R,E.D 


=  Si  • 


(4.2.22) 


Ri  +  Si  =  1  . 


(4.2.23) 


Ri  1  +  EiD  i_  )  “  1  . 


(4.2.24) 


i  b tt  * 

1  +  E.D  1 


(4.2.25) 


Si  - 


1  +  E.D 


(4.2.26) 


Finally 


Fi<V  ■  — ' 'Eb  .!  S'  I+DVl)  W 

1  +  E.D  1  ' 


(4.2.27) 


On  the  basis  of  fallout  observations,  D  has  been  taken  as  100  um. 
(Subroutines  FRATIO  and  GXPSR) 
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4.3  APPROXIMATE  ACTIVITY  CALCULATION 


For  a  pure  air  burst  (A  >  180  ft  KT”1/3*4)  or  for  a  particle  distribu¬ 
tion  that  is  specified  to  be  a  size-activity  distribution,  the  code  performs 
a  conventional  calculation  in  terms  of  a  K  factor  (Roentgens  m2  KT" 1  hr-1) 
and  a  t"1*26  decay  factor  (t  in  hours). 

If  F  is  the  fraction  of  total  cloudborrte  activity  contained  in  a  deposit 
increment  of  fallout,  the  exposure  rate  or  exposure  at  three  feet  above  one 
square  meter  of  smooth  ground  over  which  the  fallout  is  uniformly  spread 
(Roentgens  m2  hr-1  or  Roentgens  m2)  is 

A  =  FWfK$  (4.3.1) 

where  Wp  is  total  fission  yield  (KT),  and  for  exposure  rate 

$  =  t-1*26 

while  for  exposure  from  ti  to  tz 


$  =  (ti0-26  -  t£°-26)/0.26  . 

(Subroutines  PAM1A,  PAM2A  and  PCHECK) 

4.4  HEIGHT-OF-BURST  EFFECT 

Fcr  scaled  heights  of  burst  in  the  range  0  <  A  <  180  (ft  KT-1/3*4)  there 
is  at  present  not  a  satisfactory  means  available  to  compute  the  cloudborne 


particle  size  distribution  nor  to  distribute  activity  over  this  distribution. 
Complex  calculations  would  be  required  to  develop  the  information  required 
and  this  work  has  not  been  done. 

DELFIC  uses  a  very  rough  procedure  to  correct  for  height  of  burst. 

This  makes  use  of  empirical  data  that  relate  fraction  of  total  activity  ob¬ 
served  in  local  fallout  to  scaled  height  of  burst.35  In  terms  of  the 
fraction  of  local  surface  burst  activity  observed,  fj,  Showers36  has  ex¬ 
pressed  these  data  in  the  form 

fd  «  (0.45345)X/65  .  o  <  x  (4.4.1) 

where  x  is  scaled  height  of  burst  in  terms  of  ft  KT"1/3.  The  fission  yield 
is  simply  scaled  down  by  this  factor  to  correct  for  height  of  burst.  Of 
course,  if  A  a  180  (ft  KT"1/3*4),  we  assume  that  no  local  fallout  occurs, 
and  the  pure  air  burst  mode  of  calculation  proceeds.  (Subroutine  0PM1) 


5.  MAP  PREPARATION 


5.1  GENERAL  DISCUSSION 

In  this  chapter  we  discuss  the  operations  of  the  Output  Processor 
Module  (0PM).  The  0PM  receives  fallout  deposit  increment  descriptions  from 
the  Diffusive  Transport  Module  and  processes  them  into  fallout  maps  speci¬ 
fied  by  the  user.  Activity  calculations  are  performed  as  required  for  the 
requested  maps.  The  maps  are  printed  by  a  standard  line  printer.  A  printed 
and  punched  card  output  also  can  be  prepared  w.'ich  consists  of  coordinates 
of  points  on  specified  contours  in  the  map*. 

Eighteen  map  types  (Table  1),  of  which  sixteen  are  unique,  may  be  re¬ 
quested.  All  but  three  of  these  are  available  for  any  of  the  three  levels 
of  activity  calculation  (sec.  4.1);  the  exceptions  being  options  9,  10  and 
14  which  are  available  only  via  the  completely  rigorous  calculation. 


5.2  DEPOSIT  INCREMENT  PROCESSING 

In  sec.  3.4  above  we  discuss  how  an  individual  deposit  increment  (i.e., 
a  ground  impacted  fallout  parcel)  is  described  in  terms  of  its  impact  point 
coordinates,  xp,  yp,  orientation  axis  of  the  standard  deviation  ellipse,  a, 
downwind  and  crosswind  standard  deviations,  o  ,  a  ,  along  with  other  proper- 

II  I 

ties  such  as  particle  diameter  and  mass.  (See  Figs.  2  and  4.)  These  quanti¬ 
ties  are  unique  for  each  deposit  increment  and  it  is  the  principle  task  of 
the  0PM  to  compute  and  sum  contributions  from  all  deposit  increments  at  each 
map  point. 

Consider  a  deposit  increment  with  area  integrated  activity  or  total  mass 
Q.  Then  at  a  point  x,y,  the  activity  or  areal  density  of  mass  w(x,y),  is 

“(x.y)  =  2Vo(lo~  e*p 


(X-Xp)2 
2  a  2 


(V-Yp)2 

2a  2 


(5.2.1) 
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TABLE  1 


MAP  REQUEST  OPTIONS 


Map  Option 

Code,  NREQ  Description 

0  Termination  of  request  set. 

1  Count  of  fallout  deposit  increments  that  contribute  to 
each  map  ordinate. 

2  Exposure  rate  normalized  to  H  +  1  hour  (Roentgen  hr"1). 

3  Exposure  rate  at  time  H  +  T1  hours,  accounting  for  time 
of  arrival  of  fallout.  (Roentgen  hr"1) 

4  H  +  1  hour  normalized*  exposure  rate  resulting  from 
particles  in  diameter  range  T1  to  T2  micrometers.** 
(Roentgen  hr"1). 

5  Integrated  exposure  from  H  +  T1  hours  to  infinity, 
accounting  for  time  of  arrival  of  fallout  by  the  approxi¬ 
mate  nrthod.+  (Roentgen). 

6  Integrated  exposure  from  H  +  T1  to  H  +  T2  hours,  account¬ 
ing  for  time  of  arrival  of  fallout  by  the  approximate 
method. +  (Roentgen) 

7  Integrated  exposure  from  H  +  T1  to  H  <  T2  hours  assuming 
all  fallout  has  arrived  by  H  +  T1  hours.  (Roentgen). 

8  Integrated  exposure  from  H  +  T1  hours  to  infinity  assuming 
all  fallout  has  arrived  by  H  +  T1  hours.  (Roentgen). 

9  Integrated  exposure  from  H  +  T1  hours  to  infinity,  ac¬ 
counting  for  time  of  arrival  of  fallout  by  the  exact 
method.  +  (Roentgen) 

10  Integrated  exposure  from  H  +  T1  to  H  +  T2  hours,  account¬ 
ing  for  time  of  arrival  of  fallout  by  the  exact,  method. ++ 
(Roentgen). 

11  Mass  of  fallout  per  unit  area  (kg  m“3). 


TABLE  1  (con't.) 


Map  Option 

Code,  NREQ 

Description 

12 

Mass  of  fallout  per  unit  area  deposited  from  H  +  T1  • 
H  +  T2  hours  (kg  m-3). 

to 

13 

Mass  of  fallout  per  unit  area  deposited  by  particles 
diameter  range  T1  to  T2  micrometers.**  (kg  nf 3) 

in 

14 

Activity  per  unit  area  from  an  individual  mass  chain  at 

T1  hours  in  units  of  curies  m"2,  or  in  equivalent  fissions 
m"2  if  T1  =  0. 

15 

Time  of  onset  of  fallout,  (s) 

16 

Time  of  cession  of  fallout,  (s) 

17 

Diameter  of  smallest  particle  deposited,  (pm) 

18 

Diameter  of  largest  particle  deposited,  (pm) 

A  "normalized"  calculation  is  one  in  which  it  is  assumed  that  all  fallout 
is  deposited  by  H  +  t  regardless  of  actual  deposition  time. 

'k'k 

When  specifying  T1  and  T2  particle  diameters,  make  T1  slightly  smaller  and 
T2  slightly  larger  than  the  tabulated  central  diameters  for  the  particle 
size  classes. 

The  t"1*26  decay  factor  is  used  to  compute  exposure  rate  at  times  other 
than  H  +  1  hour  (sec.  4.3),  though  activity  at  H  +  l  hour  may  be  calculated 
by  the  rigorous  method.  (See  sec.  4.1) 

++  Warning:  This  calculation  probably  will  consume  a  lot  of  computer  time. 

A  complete  activity  calculation  is  done  for  each  deposit  increment  of  fall¬ 
out.  Consider  using  one  of  the  approximate  methods  (requests  5  and  6). 
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where 


X  =  x  cosct  +  y  sina 


Y  =  y  cosa  -  x  sina 


(5.2.2) 


and  X  and  Y  are  defined  similarly.  The  x,y,  and  X,Y  coordinate  axes  are 
P  P 

related  as  shown  in  Fig.  4. 


Figure  4.  A  deposit  increment  concentration  ellipse 
in  the  x,y  plane. 
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A  parameter  wm^n»  specified  by  the  user  or  calculated  by  the  program, 

(II,  Appendix  B),  is  defined  which  represents  the  lower  threshold  value  for 

contributions  from  individual  deposit  increments;  that  is,  any  u-,(x,y)  <  u  . 

min 

is  ignored.  To  determine  the  elliptical  envelope  within  which  o»(x,y)  a  wn^n» 
we  define  a  quantity  y 


Y 


y  >  0 


From  eq.  (5.2.1)  it  is  apparent  that 


(5.2.3) 


Y 


(X-Xp)g  ( V-Yp)2 

2o  2  2o  2 


(5.2.4) 


which  is  the  equation  of  the  desired  ellipse. 

Figure  4  shows  such  an  ellipse  with  tangent  lines  drawn  parallel  with  the 
x,y  axes  and  cutting  these  axes  at  points  Xy.xj,  y j . y-j- .  From  the  general  prop 
erties  of  an  ellipse,  it  can  be  shown  that  these  lines  are 


XT  *XT 


x  ±  / 2y  (a2  cos2a  +  a2  sin2a)  , 

P  II  ^ 


(5.2.5) 


and 


y 


T 


sin2a  +  o2  cos2a)  . 

I 


(5.2.6) 


In  the  code,  computation  of  x-pX-j.  and  yj»y|  is  used  to  establish  whether  or 
not  a  deposit  increment  contributes  to  a  particular  map  or  map  section. 
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For  each  deposit  increment,  map  points  are  considered  row-by-row.  The 
bounding  rows  are  determined  from  the  yy.yj  values.  For  a  particular  row. 


the  bounding  x  coordinates,  x  ,x',  are  given  by 


(y-y  )(- —  —  ]  sina  cosa  ± 
p  \a2  a2 I  * 


x  .x'  =  X  + 
c  c  p 


C0S2a  .  sin2a 


(5.2.7) 


In  preparing  maps  for  time  of  onset  or  cessation  of  fallout,  or  largest 
or  smallest  particle,  the  procedure  above  is  applied  with  u)mi-n  for  mass  per 
unit  area  of  fallout. 


6.  VALIDATION 


6.1  DISCUSSION 

★ 

Predictions  are  compared  with  observed  H  +  1  hour  normalized  expo¬ 
sure  rate  maps  for  the  five  test  shots  described  in  Table  2.  The  pre¬ 
dictions  were  executed  as  discussed  in  reference  38,  using  the  data 
listed  there.  Observed  fallout  patterns  were  taken  from  DASA  1251. 

Three  methods  of  comparison  of  fallout  patterns  are  used  here: 

1.  Visual  comparison  of  contour  maps. 

** 

2.  Comparison  of  contour  areas,  and  hotline  lengths  and  azimuths. 

3.  The  Rowland-Thompson  Figure-of-Merit  (FM).39  (Appendix  C) 

These  are  roughly  in  order  of  importance. 

Statistical  data  are  in  Table  3  and  the  contour  plots  are  on  pp.  64 
through  73.  The  contours  were  drawn  by  a  30- inch  Cal  comp  plotter,  and 
each  observed-predicted  pair  are  to  the  same  scale. 

’ *  •  I, 

TABLE  2 

TEST  SHOT  DATA 


Shot 

Total 

Yield 

(KT) 

Fission 

Yield 

(KT) 

HOB 

(m) 

Altitude 
of  GZ 
(m) 

Site 

Johnie  Boy 

0.5 

0.5 

-0.584 

1570.6 

NTS+ 

Jangle-S 

1.2 

1.2 

1.067 

1284.7 

NTS 

Small  Boy 

low 

- 

3.048 

938.2 

NTS 

Koon 

150. 

- 

4.145 

0.0 

Bikini 

Zuni 

3380. 

- 

2.743 

0.0 

Bikini 

+Nevada  Test  Site 

A  "normalized"  exposure  rate  map  is  constructed  on  the  assumption  that  all  local 
fallout  is  down  at  the  specified  time,  regardless  of  its  actual  deposition  time. 

**Hotline  length  is  defined  as  the  furthest  distance  from  ground  zero  on  a  contour, 
and  hotline  azimuth  is  the  angle,  measured  clockwise  from  north,  to  the  point  of 
furthest  distance  from  ground  zero  on  a  contour. 
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TABLE  3 

COMPARISON  OF  OBSERVED  AND  PREDICTED  FALLOUT  PATTERN  STATISTICS 

Observed/Predicted 


Test  Shot 

FM 

Contour 

(Roentgen  hr”1) 

Area 

(km2) 

Hot) ine 
Lenqth(km) 

Azimuth 

(deg) 

Johnie  Boy 

0.182 

1000 

0.278/0.029 

1.38/0.32 

359/0 

100 

0.539/0.774 

2.73/2.58 

345/344 

50 

1.271/1.787 

4.10/4.13 

343/343 

58(42)* 

28(3)* 

Jangle-S 

0.483 

500 

0.117/0.144 

0.59/1.00 

342/353 

300 

0.386/0.316 

1.50/1.23 

346/354 

100 

1.437/2.242 

3.74/5.87 

1/355 

35 

3.114/5.077 

5.06/7.68 

6/355 

40(45) 

43(42) 

Small  Boy 

0.308 

1000 

0.216/0.047 

1.00/0.25 

71/66 

500 

0.528/0.135 

1.62/0.56 

73/80 

200 

0.942/0.564 

2.22/1.69 

72/73 

100 

3.75/1.10 

5.66/3.72 

72/74 

50 

9.03/4.38 

8.10/6.47 

75/72 

63(59) 

44(36) 

Koon 

0.287 

500 

32.0/26.0 

10.2  / 12 . 5 

18/0 

250 

122/87.3 

17.3  / 24 . 2 

15/4 

100 

550/261 

41.0  /39. 5 

17/3 

33(40) 

22(22) 

Zuni 

0.105 

150 

474/2239 

98/78 

12/337 

100 

2761/3619 

125/96 

17/337 

50 

6187/6660 

138/121 

27/338 

30 

10950/9913 

177/153 

33/340 

105(16) 

17(16) 

Mean  absolute  percent  errors.  The  value  in  parentheses  is  calculated 
without  including  the  data  for  the  highest  activity  level  contour.  See 
footnote  next  page. 
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Prediction  accuracy  is  seen  to  be  good,  particularly  for  the  low  yield 

* 

shots.  Overall  mean  absolute  percent  errors  for  contour  area  and  hotline 
length  are  61  and  32  percent  respectively.  If  the  data  for  the  highest 
activity  level  contours  are  excluded,  we  have  42  and  26  percent  for  these 
errors.  The  highest  level  contours  are  particularly  difficult  to  predict, 
usually  being  in  the  region  affected  by  throwout  and  induced  activity  in 
and  around  the  crater.  DELFIC  does  not  address  this  portion  of  the  activity 
field  since  fallout  is  a  negligible  contributor  to  casual ities  there. 

It  is  important  to  emphasize  that  this  level  of  prediction  competency 
is  achieved  without  a  posteriori  adjustment  or  calibration  of  any  aspect  of 
the  model  such  as  to  improve  agreement  with  any  observed  fallout  pattern. 

The  three  low  yield  shots  were  executed  at  the  Nevada  Test  Site,  and 
their  fallout  patterns  were  measured  over  land.  For  this  reason,  observed 
patterns  for  these  shots,  though  not  highly  accurate,  may  be  considered  to 
be  superior  to  the  patterns  of  the  high  yield  shots  which  were  executed  on 
Bikini  Atoll  in  the  South  Pacific.  Not  only  are  the  fallout  fields  of  the 
high  yield  shots  very  large,  which  adds  to  measurement  problems,  but  most 
of  the  fallout  from  these  shots  fell  into  water.  Even  so,  most  of  the  Koon 
pattern  area  was  covered  by  an  array  of  fallout  collection  stations,  so  this 
pattern  is  probably  reasonably  accurate.  Zuni,  on  the  other  hand,  is  a 
special  case.  The  fallout  pattern  used  here  is  exclusively  downwind  of  the 
atoll  and  was  determined  by  an  oceanographic  survey  method  that  was  known 
to  be  inaccurate.  The  close-in  pattern  in  the  region  of  the  atoll  is 
available,  but  contains  no  closed  contours  so  could  not  be  used  here;  thus 
the  high-activity  portion  of  the  observed  pattern  for  this  shot  is  ignored, 
and  this  alone  must  account  for  a  substantial  portion  of  the  disagreement 

For  n  observed-predicted  data  pairs,  mean  absolute  percent  error  is 

n 

100  .  I  / 

n  2^\  xobs,i  "  xpred,i  I'Xobs,i 
i  =  l 
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between  observation  and  prediction  for  this  shot,  particularly  with  regard 
to  contour  areas  and  overlap  (Table  3).  In  addition,  we  have  the  following 
problem. 

Predictions  for  these  high  yield  shots  are  expected  to  be  inferior  to 
those  for  these  low  yield  shots.  This  is  because  both  of  the  high  yield 
shots  were  detonated  over  coral  soil,  and  in  the  case  of  Zuni,  a  large  but 
uncertain  amount  of  sea  water  was  lifted  by  the  cloud.  The  particle  size 
distribution  used  for  these  predictions  is  typical  of  fallout  produced 
from  the  siliceous  soil  found  at  the  Nevada  Test  Site.  We  have  not  suc¬ 
ceeded  in  developing  a  distribution  appropriate  for  coral  and  coral-sea 
water  mixtures. 

More  details  concerning  the  prediction  calculations  and  test  shot 
characteristics  are  in  reference  38. 
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APPENDIX  A 


THE  LOGNORMAL  DISTRIBUTION 


FUNDAMENTALS 

A  variable  x  is  said  to  be  normally  distributed  if  the  probability  of 
its  occurrence  in  the  range  x  to  x  +  dx  is  given  by 


dN(x|u,a2) 


(A.  1 ) 


where  u  is  the  mean  value  of  x  and  a2  is  the  variance  of  x.  The  square  root 
of  the  variance,  a,  is  called  the  standard  deviation. 

To  define  a  lognormal  distribution,  we  make  the  transformation 


x  =  i.n(y) 


(A.  2) 


In  terms  of  the  variable  y  eq.  (A. 1)  becomes 

dn(y|u,o2)  =  exp  j^-  |(^0~  u)2j  d(*ny)  ,  (A. 3) 

and  y  is  said  to  be  lognormally  distributed.20 

Some  statistical  properties  of  y  are  as  follows: 


mean{y) 

*  exp  (y  +  |  o2) 

(A. 4) 

median(y) 

«  exp(y) 

(A.5) 
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which  can  be  shown  to  be20 


n(y|y  >c2).!  =  n  (y  |  y  +  ja2,a2)  .  (A. 13) 

J 

The  moment  distributions  provide  simple  relationships  between  lognormal  dis 
tributions  of  number,  surface  area,  and  volume  of  particles  with  respect  to 
their  diameters. 


APPLICATION  TO  PARTICLE  DISTRIBUTIONS 

In  discussions  of  lognormal  particle  distributions,  confusion  frequently 
arises  because  distinction  is  not  clearly  made  between  p  and  y  and  between  a 
and  s.  Since  numerical  values  of  y  and  a  depend  on  the  base  of  the  logarithms 
used,  we  confine  our  discussion  to  y  and  s. 

Suppose  we  have  plotted  cumulative  numbers  of  particles  versus  diameter 
on  log-probability  graph  paper  and  have  obtained  the  curve  shown  in  Fig.  A.'l. 
This  straight-line  curve  indicates  that  the  distribution  of  particle  number 
with  respect  to  diameter,  <5,  is  lognormal.  Thus,  <5  is  equivalent  to  y  in 
eq.  (A. 3)  and  from  eqs.  {A. 8)  and  (A. 9)  we  have 


r  6S0 


and 


5  "  <5 bu . 1 3/ 5 bo  or  s  5  Sso/Sis.s?  • 


These  are  the  quantities  DMEAN  (urn)  and  SD  (dimensionless),  respectively,  that 
are  required  as  input  to  subroutine  I  CM.  (DMEAN  is  the  median  particle  di- 
awtee  in  th«*  distribution  o ♦  numbers  of  particles  w'th  respect  to  their 
c  1  *«w?  t  **«■  s 


PARTICLE  DIAMETER,  e(um) 


As  noted  above,  the  properties  of  the  moment  distributions  are  useful 
in  interrelating  distributions  of  particle  number,  surface  area,  and  volume 
with  respect  to  particle  diameter.  This  is  because  the  number  distribution 
is  the  zeroth  moment  distribution  with  respect  to  diameter,  surface  area  is 
distributed  via  the  second  moment  distribution,  and  volume  is  distributed  via 
the  third  moment  distribution.  Thus,  if  we  assume  spherical  particles  and  if 
the  parameters  y  and  a  are  known  for  either  the  particle  number,  or  particle 
area,  or  particle  volume  distribution  with  diameter,  then  the  other  distribu¬ 
tions  can  be  determined  from  the  equations  below.  The  parameter  a  is  the 
same  for  all  three  distributions.  If  we  use  N,  S,  and  V  as  subscripts  to 
denote  number,  surface  area,  and  volume,  respectively,  we  have  from  eq.  (A. 13) 

ys  =  +  2a2 

yv  =  uN  +  3a2 

where  y  and  a  are  related  to  y  and  s  by  eqs.  (A. 8)  and  (A. 9). 

If  base  10  logarithms  are  used  instead  of  natural  logarithms,  we  dis¬ 
tinguish  the  distribution  parameters  by  use  of  primes,  y^  and  a',  and  the 
relations  become 

v's  =  y‘  +  2zn(10)(a‘)2 

Uy  »  UjJ|  +  3£n(10)(a'  )2 

where  zn{10)  =  2.3026. 

The  distribution  of  particle  nass  with  respect  to  diameter  is  taken  to 
be  equivalent  to  the  volume  distribution  with  respect  to  diameter.  This  im¬ 
plies  that  ali  particles  have  the  same  density. 


DISTRIBUTION  HISTOGRAM 


For  computation  purposes,  the  continuous  lognormal  distribution  is 
replaced  by  a  histogram.  The  computer  program,  via  subroutine  DSTBN,  does 
this  automatically  by  use  of  the  distribution  parameters  and  the  number  of 
size  classes,  NDSTR,  which  is  input  by  the  user  or  on  default  of  input  is 
set  to  100. 

lijj,  a,  and  uv  are  determined  from  DMEAN  and  SD  via 
yN  =  nn (DMEAN) 
o  =  &n(SD) 
uv  =  pN  +  3a2  . 

Define  the  normal  distribution  function  argument  x  as 

£n(6)  -  Uy 


where  6  is  particle  diameter.  Then 

x 

N(x)  =  f  exp(-t2/2)dt 

/2v  . 

Subroutine  DSTBN  constructs  the  particle  size  class  histogram  as  follows. 
Each  size  class  contains  a  constant  volume  fraction,  ANy,  of 


AN, 


H2 


1/NDSTR 


Let  D.j ,  i  =  1,2 . ,NDSTR,  be  the  upper  (i.e.,  the  larger  particle) 

boundary  diameter  of  the  i-th  particle  size  class.  The  histogram  is 
ordered  with  the  largest  particles  in  the  first  size  class,  and  so  on. 
Then,  for  the  i-th  size  class 

N(Xi)  =  I'ANy 


and 


an(D.+1)  =  x.a  +  . 

The  upper  boundary  of  the  first  size  class,  Dl9  and  the  lower  boundary 
of  the  last  size  class,  dndstr+1 *  are  sPecia^  cas&s.  These  are  taken  to  be 
the  diameters  at  ANy/2  and  l-ANy/2,  respectively.  That  is, 


and 

ANy 

N^XN0STRvl^  =  1  '  T 


In  these  calculations  x  is  determined  from  given  N(x)  via  equation 
26.2.23  of  ref.  37. 

The  central  particle  diameter  for  the  i-th  class,  6..,  is  given  by 


<5 


i 


/dTd 


i+i 


H  3 


If  NDSTR  =  1,  a  single  size  class  is  created  with 
Di  =  (DMEAN)  *  (5.0  *  SO) 

02  =  (DMEAN)  /  (5.0  *  SD) 


«l  = 


DMEAN  . 


APPENDIX  B 

POWER  LAW  PARTICLE  SIZE  DISTRIBUTION 


DATA  ANALYSIS 

Suppose  that  we  have  a  sample  of  fallout  particles  of  mass  m  and  average 
density  p,  and  we  size  the  sample  into  N  fractions,  the  ith  fraction  of  mass 
AF.j  containing  particles  in  the  diameter  range  A6^  centered  on  <5^.  To  obtain 
the  power  law  distribution  parameters  k/m  and  X  (sec.  2.1.6),  we  plot 
loga (aF./aa . )  vs.  log  (6.)  where  a  is  the  logarithm  base.  A  straight  line 
is  fitted  to  the  plot,  and  we  obtain  for  this  line 


intercept  =  c  =  "log.^^5 


slope  =  s  =  3  -  X  . 


Therefore,  the  DELFIC  distribution  parameters  are 


X  =  3  -  s 


k/m  =  6  ac/(7rp) 


where  m  is  in  kg,  <5  in  meters  and  p  in  kg  m"3. 


DISTRIBUTION  HISTOGRAM 

DELFIC  creates  a  histogram  representation  of  the  power  law  distribution 
(Subroutine  DISTBN)  for  use  in  the  fallout  calculations.  The  histogram  con¬ 
sists  of  NDSTR  particle  size  classes  of  equal  mass  fraction  AF  -  1/NDSTR. 


Let  Di  be  the  upper  (i.e.t  larger  particle)  boundary  of  the  i-th 
particle  size  class.  The  table  is  ordered  with  the  largest  particles  in 
the  first  class,  and  so  on.  If  we  assume  that  the  smallest  particle  in 
the  NDSTRth  class  is  much  smaller  than  DNDSTR,  we  see  from  Eq.  (2.1.18) 
that 


0 


4-X 

NDSTR 


6(4  -  X)AF  /nu 
■np  'k' 


By  recursive  use  of  this  relation  with  Eq.  (2.1.18),  we  find  that 

1 

0,  ■  (nostr  -  1  +  1)4-'X  dnostr  . 


Size  class  central  diameters,  6..,  are 


5i 


■  -TC 


To  establish  a  central  and  lower  boundary  diameter  for  the  NDSTRth  class, 
we  say  that 


and 


6NDSTR 


dndstr+i 


^ndstr^2/dndstr  • 
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APPENDIX  C 


FALLOUT  PATTERN  COMPARISON  BY  THE  F I GURE-OF-MER IT  METHOD 


Rowland  and  Thompson39  developed  this  method  for  comparison  of  pairs 
of  fallout  contour  maps  by  computation  of  a  single  index,  the  FM,  that  is 
a  measure  of  contour  overlap  between  them.  For  each  contour  common  to  the 
patterns,  the  area  overlapped  and  the  area  not  overlapped  is  calculated. 

The  areas  are  weighted  by  the  average  radiation  level  between  successive 
contours.  Sums  over  all  contours  of  weighted  overlapped  areas  and  weighted 
total  areas  are  computed,  and  the  FM  is  the  ratio  of  the  two  sums.  For  com¬ 
pletely  overlapped,  perfectly  matched  patterns,  FM  =  1;  for  no  overlap,  FM  - 
Mathematically,  FM  is 


h  -  -M> 

i 'Wi) 


where 

N  is  the  number  of  contours  in  the  patterns.  The  summations 
are  from  highest  contour  to  lowest 


r.  is  activity  of  the  i  th  contour  (Roentgen/hr),  rQ  =  10  rj 

a..  is  common  (i.e.,  overlapped)  area  for  the  i  th  contours. 

a  =  0. 
o 

A.,  is  total  area  of  the  i  th  contour.  The  summation  in  the  de¬ 
nominator  is  computed  for  both  patterns,  and  the  largest  sum 
is  used.  A  =  0. 
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The  FM  has  been  found  to  have  limited  utility  as  a  measure  of  fallout 
prediction  accuracy.  This  is  mainly  for  two  reasons.  First,  and  most 
important,  is  that  being  a  measure  of  overlap,  the  FM  is  stronqly  biased 
in  favor  of  overprediction;  that  is,  it  favors  predictions  that  cover  a 
large  area,  and  therefore  overlap  the  observed  pattern,  regardless  of 
other  considerations.  Second,  the  FM  method  imposes  no  penality  for  miss¬ 
ing  or  extra  contours;  contours  not  common  to  both  patterns  are  simply 
ignored. 
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APPENDIX  D 


GLOSSARY  OF  SYMBOLS 


activity 

specific  heat  of  air  and  water  vapor  at  constant  pressure 
specific  heat  of  dry  air  at  constant  pressure 
specific  heat  of  water  vapor  at  constant  pressure 
specific  ',ieat  of  soil 

scaled  depth  of  burst  (ft  KT"1/3*4), differential  operator 
particle  diameter  at  size  class  boundary 
turbulent  kinetic  energy  per  unit  mass  of  cloud 


particle  settling  speed 


rate  of  fallout  from  the  cloud 


total  number  of  equivalent  fissions 
acceleration  of  gravity 


heat  energy 


vertical  cloud  radius 


relative  humidity 

Von  Karmai.'-i  constant  (0.35) 
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Glossary  (con't.) 


k2 


k3 

L 


m 


m. 


m. 


m.. 


Nr 


ws 


P 

pv 

q(x) 

Q 


cloud  rise  kinetic  energy  to  turbulent  energy  conversion 
parameter  (eq.  (2.2.1a)) 

cloud  turbulent  energy  dissipation  constant  (eq.  (2.2.4a)) 

latent  heat  of  vaporization  of  water  or  ice,  Monin- 
Obukhov  length 

mass  of  cloud 

mass  of  fallout 

mass  of  air 

mass  of  water 

Davies  number  (eq.  (2.2.9)) 
atmospheric  pressure 
saturation  water  vapor  pressure 

,  where  x  is  water  vapor  mixing  ratio 

A  ■  X 

activity  or  mass  in  a  fallout  parcel 
horizontal  position  vector 
ideal  gas  law  constant  for  air 
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Glossary  (con't. ) 


cloud  radius 

geometric  standard  deviation  of  particle  diameter  for  a  log¬ 
normal  particle  size  distribution,  soil  mixing  ratio  (soil  mass 
per  unit  mass  of  dry  air),  particle  settling  speed  slip  factor 
(eq.  (2.2.11)),  measurement  height  of  surface  wind  speed 

4ttR2 

time 

time  of  the  second  fireball  temperature  maximum 
temperature 

Tq(x),  i.e.,  virtual  temperature 

particle  deposition  time  from  the  base  of  the  jth  wind  stratum 
cloud  rise  velocity 

air  friction  velocity  in  the  atmospheric  surface  layer  (eq. 
(3.3.4)) 

horizontal  wind  vector 

characteristic  cloud  rise  velocity  (eq.  (2.2.1b)) 
cloud  volume 

condensed  water  mixing  ratio  (mass  of  condensed  water  per 
unit  mass  of  dry  air),  vertical  wind  component 


Glossary  (con't.) 


total  energy  yield  of  the  nuclear  explosion  (KT) 
fission  yield  (KT) 

water  vapor  mixing  ratio  (water  vapor  mass  per  unit  dry  air 
mass),  coordinate  in  the  west-to-east  direction 

power-law  particle  size  distribution  exponential  parameter 

fission  yield  of  jth  mass  chain 

vertical  coordinate 

aerodynamic  surface  roughness  length  (eq.  (3.3.4)) 

wind  direction  angle,  limiting  vertical  distance  for  wind 
data  interpolation  (eq.  (3.5.2)) 

limiting  horizontal  distance  for  wind  data  interpolation 
(cq.  (3.5.2) ),  ratio  of  Lagrangian  to  Eulerian  time  scales 
(eq.  (3.3.2)) 

l+x 

ratio  of  gas  density  to  total  density  of  cloud  =  y+Yfs+w 

particle  diameter 

median  particle  diameter 

turbulent  energy  density  dissipation  rate 


air  viscosity 


Glossary  (con't.) 


X 

A 

y 

5 

p 

p 

a 


w(x,y) 


0) 


min 


SUBSCRIPTS 

a 

B,b 

e 


scaled  height  of  burst  (ft  KT"1/3),  activity  decay  constant 

scaled  height  of  burst  (ft  KT”1//3‘4) 

cloud  entrainment  parameter  (eq.  (2.2.5a)) 

ratio  of  molecular  weights  of  air  and  water  (29/18) 

air  density 

particle  density 

standard  deviation,  o2  variance 

wind  shear  stress  in  the  atmospheric  surface  layer  (eq.  (3.3.3)) 

fraction  of  available  heat  energy  in  the  initial  cloud  used  to 
heat  air  and  soil 

contribution  of  activity  or  mass  from  a  deposit  increment  of 
fallout  at  map  point  x,y 

threshold  activity  or  mass  contribution  from  a  single  fallout 
deposit  increment  at  a  map  point 


a  ir 

cloud  or  fallout  parcel  base 
ambient 


ert 


entra inment 


Glossary  (cont.) 


SUBSCRIPTS  (con't.) 

g,GZ  ground 

i  initial  cloud 

r  refractory  mass  chain 

s  stabilization,  soil 

T,t  cloud  or  fallout  parcel  top 

w  water  or  water  vapor 

n  along  wind 
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